Submitted to Physical Review B 



Phase Transitions in the Two-Dimensional XY Model 
with Random Phases : a Monte Carlo Study. 

J. Maucourt, and D.R. Grempel 

CEA/Departement de Recherche Fondamentale sur la Matiere Condensee. 
SPSMS, CENG, 17, rue des Martyrs, 38064, Grenoble Cedex 9, France. 

(February 1, 2008) 



Abstract 

We study the two-dimensional XY model with quenched random phases by 

Monte Carlo simulation and finite-size scaling analysis. We determine the 

phase diagram of the model and study its critical behavior as a function of 

disorder and temperature. If the strength of the randomness is less than a 

critical value, o~ c , the system has a Kosterlitz-Thouless (KT) phase transition 

from the paramagnetic phase to a state with quasi-long-range order. Our data 

suggest that the latter exists down to T = in contradiction with theories that 

predict the appearance of a low-temperature reentrant phase. At the critical 

disorder Tkt — ► and for a > a c there is no quasi-ordered phase. At zero 

temperature there is a phase transition between two different glassy states at 

a c . The functional dependence of the correlation length on a suggests that 

this transition corresponds to the disorder-driven unbinding of vortex pairs. 
PACS numbers: 75.10.Nr, 75.40.Mg, 74.40 
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I. INTRODUCTION 



The two-dimensional XY model with quenched phase disorder describes the thermody- 
namics of a variety of systems of experimental interest. Examples of these include magnetic 
systems with random Dzyaloshinkii-Moriya interactions ffl, crystal systems on disordered 
substrates 0, arrays of Josephson junctions with positional disorder and vortex glasses 
0. The critical properties of these systems are described by the model Hamiltonian 

H=-J cos(^-^-Ai), (1) 

<i,3> 

where the sum runs over the Nj, bonds of a 2 — d square lattice, 0$ is a phase variable attached 
to site i and Ay = —Aji is a quenched random bond-variable whose precise physical meaning 
depends on the nature of the system described by the model. We assume for simplicity that 
the phase shifts on different bonds are uncorrelated and that each configuration {A^} occurs 
with probability 

I 

— — exo — ■ 

(27TCT 

The nature of the phase diagram of this system in the a — T plane has been controversial. 
By mapping model ([I]) into a two-dimensional Coulomb gas in a random background of frozen 
dipoles, Rubinstein et al. [|TJ generalized the original treatment || of the Kosterlitz-Thouless 
(KT) transition to the random case. Their analysis leads to the phase diagram shown in 
Fig. [I]. For a > a c = n/8 the system is paramagnetic at all temperatures. For a < er c , in 

define the 



"27 ^ A ^ 

zo <ij> 



(2) 



the limit of infinite vortex core energy , two lines, T±(a) = 2 j- 1 ± yl — 8a /tt 
boundaries of a quasi-ordered KT phase. The transition along the upper curve is similar 
to the KT transition of the pure system except that the transition temperature and the 
value of the jump of the exponent rj at the transition are reduced by the randomness. The 
second transition line has no counterpart in the pure case and it signals the reentrance of 
the disordered phase at low temperature. It occurs because of the presence of a term in the 
renormalizat ion-group flow equations M that makes the vortex fugacity grow continuously 



2 



with increasing length scale at low temperature. The effect of a finite vortex-core energy is 
to push the two transitions down to lower temperatures, T KT (a) < T + (a), T r (p) < T_(ct). 

So far, the reentrant phase has not been seen either in numerical simulations [^,|J or 
experiments |9j on disordered arrays of Josephson junctions. Although finite-size ]7J or 
pinning [H effects can account for the failure to observe the reentrance transition, it has been 



recently suggested flH||llfl that, in fact, the latter does not exist. The physical argument 
supporting this claim is that topological defects can appear at T = only if the energy 
cost to create them is balanced by an energy gain coming from the random field. However, 
a simple probabilistic estimate [|IIJ] shows that, in an infinite system, the probability of 
finding sites on which it pays to create an isolated vortex vanishes below a c . On the basis of 
this reasoning one expects the quasi-ordered phase to be stable below the critical disorder 
at T = 0. Refining this heuristic considerations by taking into account the interactions 



between defects, Natterman et al. |TT[ have derived a set of renormalization group equations 
that reduce to those derived earlier [jXJ] above a crossover temperature T* « 2Ja, but are 
different from them below it. In the modified theory, vortices are irrelevant in the whole 
region a < a c , T < Tkt{&), where the correlation functions exhibit power-law decay. At 
T = they find a disordered-driven transition at o = o c at which there is a jump of the 
correlation function exponent from a finite value to zero with Arj c = 1/16. Above a c the 
correlation length is finite. 

In this paper we present the results of extensive Monte Carlo simulations of the XY 
model with random phases. The main difference between our simulations and those per- 
formed previously |7|,|8| resides in the use of finite-size scaling in the analysis of the results. 
Considerations of scale invariance and universality provide a sensitive tool to distinguish 
and characterize the different possible phases by imposing stringent conditions on the form 
of the correlation functions as functions temperature, randomness, and system-size. 

Scaling behaviour in the paramagnetic phase, near the KT transition, is most con- 
veniently discussed in terms of the Binder function [|16||, essentially a ratio of moments 



of the order parameter. Parametrizing the correlation length in the KT form, £(T) 



exp Aj yT/TxT — 1, and adjusting A and Tkt so that scaling holds, we determined the 
a-dependence of the KT transition temperature and of the expected [II] reentrant transition 
temperature. We found that Tkt{°~) vanishes sharply at a c = tt/8, in agreement with a 
prediction of Ozeki and Nishimori JT2|. The measured disorder-dependence of Tkt near 
the critical disorder is consistent with the expression that can be derived from the modi- 



fied renormalization group equations [11|. We found no trace of the expected signature (cf. 
Section II A ) of a reentrant transition in the size dependence of the Binder function at low 
temperatures. 

For T < Tkt, analysis of the size-dependence of the moments of the order parameter 
gives access to the temperature dependence of the correlation function exponent, r](T,a). 
It can be shown that, in the case of a reentrant transition, drj/dT must change sign going 
from positive near Tkt to negative near T r . We found instead that i](T) is monotonic and 
assymptotically approaches the spinwave result |TTJ r](T) ~ 4~(T/ J+cr) at low temperature. 
This suggests that the quasi-ordered phase is stable down to T = for sufficiently weak 
disorder. 

A scaling analysis of the type described above is not possible above a c because in this 
region the disorder and temperature dependence of the correlation length is not known a 
priori. However, from the analysis of the data in the weak-disorder region one can determine 
not only the critical temperature and exponents but also the full shape of the scaling function. 
Since by universality the latter is the same throughout the paramagnetic region , we can 
use it to analize the resuls for a > o~ c as well. We found that in the strong disorder region 
the spin-spin correlation function decays exponentially at all temperatures. The correlation 
length decreases with increasing T and stays finite as T — > 0. The extrapolation of our 
results to T = is consistent with the form log£(T = 0, a) ~ (er — o" c ) _1 expected near a 
disorder- driven vortex-unbinding transition JTTJ . 
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II. METHOD 



In this Section we describe the technique and the method of data analysis that we have 
used in our simulations. 

A. Scaling Analysis. 

The analysis is based on the finite-size scaling properties of the moments of the order 
parameter, 

(3) 

d 

where N = L 2 is the total number of spins in the system and L is its linear dimension. 
Here the symbol (...) indicates the usual thermodynamic average with the Gibbs measure, 
P ~ exp(— j3H), and [. . .] d is the average over the phase shift distribution given in Eq. 0. 

In a finite system all configurations related by a global spin rotation enter in the ther- 
modynamic average with equal weight. Since the average over the global rotation angle of 
all odd moments vanishes, the first non-trivial moments are and q( 4 \ If L 3> a (a is the 
lattice spacing), we may write 

qW(L,T,a) = L-*f s d 2 rC 2 (r), 

(4) 

q^(L,T,a) = L~ 6 f s d 2 n d 2 r 2 d 2 r 3 C A (n,r 2 ,r 3 ), 

where the integrals are over the surface of the system and C n , n = 2,4, are the two and 
four-point correlation functions respectively. These are given by the rotationally invariant 
expressions 

C 2 (r) = | [<cos[0(TO-0(O)])] d) 

(5) 

C 4 (f u f 2 ,f 3 ) = | [<cos[0(ri) - 9{r 2 ) + 6(r 3 ) - 9(0)])] d . 

The scaling theory of phase transitions may be used to determine the general form of the 
moments as functions of temperature, randomness and system size. In the paramagnetic 
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phase, sufficiently close to Tkt for the system to be in the scaling regime, the two-point 
function scales as [IBj 

C 2 (f)~r-»C 2 (r/0, (6) 

where 77 is the correlation function exponent, £ = £(er, T) is the correlation length, and the 
function 62(2;) is universal. Using Eq. |6] and the corresponding expression for the four-point 
function, and may be written in the form : 

q^(L,T,a) = L-^ 
g( 4 )(L,T,a) = L-^Qi 4) (L/0, 

valid in the disordered phase. The functions and in the expressions above are 
universal functions of their argument. 

The correlation length diverges at Tkt arid stays infinite throughout the KT phase. In 
this phase the long-wavelength behavior of the system is described by renormalized spin- 
wave theory. The two and four-point functions may be easily calculated within this theory 
resulting in the expressions : 

C 2 (r) ~ r~\ 



C 4 (fi,r2,f 3 ) 



T2 \ri-rz\ 



ri r-j, ri-T2 \r2-r-j,\ 



Substituting Eq. [8] into Eq. |] we obtain the scaling form of the moments in the quasi- 
ordered KT phase : 

q W(r,T,a)~L-«QW(r]), 
q^(r,T,a)^L^Q^(r]), 

where the universal functions and given by 

qL 2) (t?) ~ Jd 2 x\x\-^, 



Q- (rj) = fd 2 Xid 2 x 2 d 2 x 3 



'J |T i-^3l 



.xi x s \xi-x 2 \ \x 2 -x 3 \_ 

depend on temperature and disorder only through the correlation function exponent 77. 
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In the first step of the analysis the unknown factor L~ v can be eliminated from Eqs. and 
|9] by working with the ratio u(L, T) = /[q^] 2 . It is customary to normalize this quantity 
so that it vanishes at high temperature and it goes to unity for a homogeneous ferromagnet 



at T = 0. The normalized quantity, known as the Binder function [IB|, is defined by 



g{L,T) = ^(3-u{L,T)), (11) 
It follows from Eqs. [F| and |9] that the Binder function obeys the scaling law 

g(L,T,a) = { ~ (12) 

[G4 V ) ,T<T KTj 

where G± are the universal scaling functions in the paramagnetic and quasi-ordered phases, 
respectively. Our discussion in the Section below is based on the following consequences of 
Eq.|TJ: 

1. In the paramagnetic phase g(L,T) depends on both temperature and system size. 
However, by choosing £ (T, a) appropriately one can make all the data collapse on a 
universal function, G + (L/£). 

2. At the KT transition temperature and in the KT phase £ is infinite. Therefore the 
curves for different sizes must merge at Tkt- In the KT phase the Binder function 
depends on temperature and disorder through a universal function of the correlation 
function exponent, G-(rj). 

3. A logarithmic plot of q^ 2k \L,T) as a function of L for different temperatures in the 
KT phase must give a series of straight lines whose common slope is —krj. 

4. If a reentrant transition occurs, the curves g(L,T) for the different sizes must split 
again at T r and remain distinct down to T — > 0. 

In addition, one can show (see Section [111 A 2| ) that, if there is a reentrant transition, the 
slope drj/dT must change sign at some temperature between Tkt and T r . This gives yet 
another signature of reentrance. 



B. Numerical method 



We have determined by Monte Carlo simulation the 2— nd and 4— th moments of the order 
parameter for systems of planar spins on an L x L square lattice with periodic boundary 
conditions. In studies of disordered systems it is particularly crucial to check that thermal 
equilibrium has been achieved before making the Monte Carlo measurements. We have 



done this using a procedure first introduced by Bhatt et Young [14] for Ising glasses and 



generalized by Ray and Moore |15j to the case of XY spin glasses. The method is based on 
the comparison between two quantities. One is the mean-square-averaged overlap (MSAO) 
of two time-delayed configurations of the same evolving system. The other one is the MSAO 
between the instantaneous configurations of two identical copies (or replicas) of the system 
that have evolved independently starting from arbitrary initial conditions. The two overlaps 



| I5| are defined by 
and 



N 



J2 cos [9i(t + t w ) -dfa 



(13) 



N 



Y>os et ] (t w ) - o?\t w ) 



(&), 



(14) 



respectively. Here, t w is a waiting time during which the systems considered here have 
evolved from the initial conditions, and t and (a) and (b) denote the time delay and the 
replicas referred to above. Let t eq and t r be the equilibration time and the equilibrium 
relaxation time, respectively. It follows from general considerations of ergodicity that 0(t + 
t w ,t w ) — * O r {t w ) when t w > t eq and t > t r . 

In Monte Carlo simulations thermal averages are replaced by time averages. Thus, to 
compute the time-delayed overlap the system is simulated for t w Monte Carlo time-steps per 
spin (MCS). The final configuration is stored, and the system is left to evolve further during 
t MCS. Measurements are taken at times t m — t + t w + m, m — 0, . . . , M, and the required 
overlap is computed as the average of the overlaps between the configurations at times t m 
and t w . Since, in general, t eq 3> t r , we may take t = t w and write |l5j : 



S 



(15) 

d 

Similarly, in order to compute O r (t w ), we simulate in parallel two copies of the system for 
t w MCS steps and we take the average of the overlaps of the next M configurations of the 
two replicas, 

(16) 

d 

In this work we have studied systems with L = 4,6, 8, 10, 12 and 16 for sixteen values of 
er in the range < a < 1 and temperatures in the range 0.3 <T/J< 1.7. For certain values 
of the disorder and for sizes L < 10 we pushed the simulations down to T/J = 0.1. The 
thermalization times vary between t w = 2 x 10 3 for the smaller systems up to t w = 2 x 10 5 
for the larger ones at low T. In all our simulations M = t w . The number of realizations of 
the random bonds simulated to perform the configuration average varies from 256 for the 
16 x 16 systems up to 2048 for the smallest ones. However, for values of o ~ a c we have 
averaged over four times us much bond configurations. Our simulations were performed on 
a 128— processor CRAY T3D parallel computer. 

III. RESULTS 
A. Weak disorder 

We show in Figs. to |2|d the numerical values of the Binder function, Eq. [11], as a 
function of temperature and system-size for four values of a, representative of our results 
below cr c « 0.393. The four curves are qualitatively similar except that the temperature 
scale shifts to the left with increasing disorder. For each size we can identify two regimes. 
At high temperatures g(L, T) decreases with increasing temperature or system-size. There is 
an inflection point (not always seen in the temperature range shown in the figures) and, for 
sufficiently high T, g(L,T) reaches a temperature-independent plateau whose hight depends 
on the system-size. At low temperatures g(L,T) is a convex function of T. Below a certain 
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temperature the Binder function is L-independent within the statistical error. There is no 
further splitting of the curves at low T in the temperature range covered by our simulations. 

These data follow closely the scenario anticipated in the paragraph below Eq. [12] for a 
system that goes from a disordered to a quasi-ordered state through a Kosterlitz-Thouless 
transition at a temperature Tkt- The latter is identified as that at which the curves for 
different sizes merge. The observed absence of size- dependence of g at low temperatures 
indicates that the domain of stability of the quasi-ordered phase extends, at least, down to 
the lowest temperatures that we simulated. It will be seen below that, except for the lowest 
concentrations, these are below the expected [Q] reentrance temperature. The quantitative 
analysis of the data is as follows. 

1. The paramagnetic region 

To discuss quantitatively the data in the high temperature region we assume that for 
each value of a the correlation length is of the KT form, 

f (T) ~ exp A/ \Jt/Tkt — 1) (17) 

where A and T/T KT are disorder- dependent constants. These are determined such that for 
each a data for all temperatures and sizes collapse into a unique function of L/£ as required 
by Eq. |12|. Examples of the resulting scaling plots for g(L, T) are shown in Fig. |3] for the same 
values of a as in Fig. [2|. It may be seen from the figure that the Binder function is a decreasing 
function of It decays exponentially for L » £ and it varies lineraly for L/£ ~ 0. At 
the transition temperature g takes the universal value g(L, Tkt) — 0.972 ± 0.001. It will be 
seen below that the differences between the curves corresponding to the different values of 
a can be absorbed in a disorder-dependent amplitude in the correlation length, Eq. [T7|. 

The dependence of the critical temperature on randomness obtained from the scaling 
analysis is shown in Fig. For weak disorder Tkt decreases slowly as a function of o. 
When the transition temperature goes below T* its variation with a becomes steeper and 
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it is seems to be very abrupt for a ~ cr c . In particular, whereas the system has a relatively 
high transition temperature for a = 0.392, the analysis of the shape of the Binder function 
(see below) shows that the T = correlation length is finite for o = 0.393. The observed 
sharpness of the transition line agrees with the prediction of Ozeki and Nishimori jT2| that, 
if a low temperature KT phase exists in this model, the phase boundary must be parallel 
to the temperature axis as T — ■> 0. It may be shown |~3[ that the modified RG equations of 
Natterman et al. fll| imply that, for o « er c , Tkt{c) ~ 2E C / In [7r 3 / (er c — a)] where E c is the 



vortex-core energy. We see from Fig. |] that our results are consistent with this expression. 

As a comparison of Figs. [| and |l] shows, Tkt is much reduced from its upper limit T + . 
This is an indication that the vortex-core energy is relatively small and may considerably 
renormalize the expected value of the reentrant transition temperature. The core energies 
can be computed from the measured Tkt by integrating the renormalization group equations 
of Rubinstein et al. along the critical trajectory. The RG equations are : 



*f = -47i 3 K 2 y 2 



(18) 

— (0 — tvK -L^K 2 rr\^, 

dl 



& ~ (2 - irK + nK 2 a)y 



where K — J/T is the running coupling constant and y is the vortex fugacity whose bare 
value is related to the temperature and E c by yo(T) = exp(—E c /T). The renormalization 
group trajectories, obtained by integration of Eq. [T3|, are given by the family of curves 

y 2 = ^-r(2/K + 7TlnK -TiaK) + C, (19) 

where C is an integration constant. The critical trajectory is determined by the condition 
that points on it flow to the fixed point y — 0, K — J/T + with T + = ^ 1 + — Sa/n 



The physical initial condition intersects the critical trajectory at two temperatures, Tkt and 
T r , solutions of : 

1 



exp(-2E c /T) 



2{K' 1 - K~ l ) + Tr In — - Tia{K - 



(20) 



2tt 3 

The vortex-core energy E c may be determined by substituting the measured values of the 
KT transition temperature in Eq. |2(| We found that E c is a decreasing function of disorder 
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that varies between 2.28J for a = and 1.9J for a c . Once the vortex-core energy is known 
the reentrance temperature may easily be computed by searching for the second solution of 
Eq. |20]. The result is shown by the open circles in Fig. |j. It may be seen that, for a > 0.25, 
T r lies above the lowest measured temperatures. The fact that the Binder function shows 
no measurable size-dependence below T r for those values of a for which the region T < T r 
is accessible to us suggests that for sufficiently weak disorder the KT phase is stable at 
low temperature. The measurements of the disorder and temperature dependence of the 
exponent r\ that we present below give further support to this conclusion. 



2. The quasi- ordered phase 

The results in the KT phase below T^t are most conveniently analyzed in terms of the 
scaling properties of the moments of the order parameter, Eqs. [5] and [10]. Although all even 
moments contain the same information, we have worked with because its temperature 
dependence is particularly simple. This is shown in Fig. |5] where we plot raw data obtained 
for a = 0.25. We see that is a decreasing function of temperature and system size that 
can be quite accurately fitted by a straight line in the temperature range of the simulations. 
According to Eq. || we expect that all the L-dependence of q^ be contained in the prefac- 
tor, L~ n . The unknown //-dependence in the function Q ( l\t]) may be eliminated from the 
problem by making plots of logg^ 4 -* vs. logL for different values of the temperature and 
of the disorder. For a fixed value of a we obtain a series of straight lines, one for each 
temnperature, whose slopes are —2r]. This is illustrated in Fig. |6|, where we represent the 
data for a = 0.25 as a function of L. 

We determined by this procedure the T-dependence of rj for all values of a. Results 
are plotted in Fig. |7| for a few values of the disorder parameter. The correlation function 
exponent increases continuously with temperature. The extrapolation of our results down 
to T = is consistent with the spinwave result, rj(T,a) = (T/J+ er)/(27r). The correlation 
function exponent at the transition, r](T KT ), is disorder-dependent and smaller than 1/4, 
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the value for the pure system. The two sets of RG equations p]Jll| give different predictions 
for r\ c . However, up to a = 0.39 the difference between the two theories is within our error 
bars and we can not decide in favour of one or the other. 

The curves shown in Fig. [7| are inconsistent with the existence of the reentrant phase. 
Since all points on the critical trajectory flow to the same fixed point, the value of the 
correlation function exponent at criticality must be the same at T r and at Tkt- This 
implies that, if there is a reentrant transition, dt]/dT must change sign at some temperature 
intermediate between Tkt and T r . Our data in the region a > 0.25 where temperatures 
below T r are accessible do not show this behaviour. 



B. Strong disorder. 

The Binder function for four values of a > a c is shown in Figs. |8|a to |]d. The behavior 
of g(L,T) in these cases is very different from that shown in Fig. [| for weak disorder. The 
Binder functions for different sizes do not collapse but they stay distinct down to the lowest 
temperatures considered. Although it is obviously impossible to guarantee that they do not 
join at lower temperatures, this seems unlikely because, in order to merge, the curves would 
have to turn upwards at low temperatures. However, in all the cases in which we did observe 
a KT transition, g(L,T) was found to be a convex function of T at low temperature. We 
interpret this behavior as evidence that for a > a c the correlation length stays finite at all 
temperatures. 

For a fixed size, g(L,T) saturates to a constant value with decreasing T indicating that, 
at low temperature, the correlation length only depends on the strength of the disorder. 
This dependence is very strong as shown in Fig. ^] where we represent the temperature and 
a dependence of g(L,T) for L = 16, our largest size, and T/J > 0.3. There is a five-fold 
decrease in the assymptotic value g(L,T) when a varies between 0.4 and 1. 

In order to analyze our results quantitatively it is convenient to fit them to a simple 
analytic expression. We have chosen a four-parameter function, 
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g(L, T) = mi + m 2 [1 - tanh[m 3 (T - T )}} , (21) 

that, despite its simplicity, accurately describes all the data throughout the paramagnetic 
phase as the solid lines in Figs. |8| and ^| show. Thanks to this parameterizatio our results 
can easily be extrapolated down to T = 0. 

The scaling analysis of the data in Figs. to |^d can not be performed in the same way 
as for a < a c because we do not have here an explicit expression giving the temperature de- 
pendence of the correlation length. However, this is not actually necessary. By universality, 
the scaling function G + (x) of Eq. [12] is the same throughout the disordered phase. If we 
determine it from the data in the weak-disorder region, we can use it in the region a > a c 
too. 

In order to do this we reconsider the scaling plots of Fig. |3|. The reason why the curves 
for different values of a are all different is that in the scaling analysis the correlation length 
has been determined up to an amplitude that depends on the strength of the disorder. We 
can compute all but one of the amplitudes by rescaling horizontally the curves of Fig. |3| 
such that they can all be superimposed to one of them. Choosing the curve for a = 
as reference and its amplitude as the unit of lenght, this procedure leads to the universal 
scaling function shown in Fig. [K| This plot where there more than six hundred points 
corresponding to different system sizes, temperatures and values of a lie on the same master 
curve is a beautiful example of universality. 

Having thus determined the scaling function G + (x)y, the correlation length in the strong- 
disorder regime can be computed by reading from Fig. [10] the values of L/£ corresponding 
to the measured values of g(L,T). This yields the temperature and disorder- dependent 
correlation length except for an undetermined scale factor. The results of the analysis are 



presented in Fig. [11]. At high temperature the correlation length is small and varies slowly 
with disorder. At low temperatures £ is T-independent and reaches a saturation value that 
increases very fast as we approach the critical disorder from the right. As seen in the figure, 
the error bars become very large when o — > a c , the uncertainties in G + (cf. Fig. |It]) and 
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in the value of the Binder function leading to large errors in £ for L/£ << 1. The results 
can be reliably extrapolated to T = only for a > 0.49. The disorder dependence of the 
ground-state correlation length in this region of values of a is represented in Fig. [12]. It can 
be seen that the data are well described by the expression £(T = 0, a) ~ exp [b/ (1 — <r c /a)] 



that follows from the T = form of the modified RG equations of Natterman et al. |TT 



This result supports their conclusion that at zero-temperature there is a disorder-driven 
transition at a c between two spin-glass states that differ by the properties of the spin-spin 
correlation function. 

In summary, we have studied the two-dimensional XY model with random phases by 
Monte Carlo simulation. An essential element in our work is the use of finite-size scaling in 
the analysis of the results. The scaling properties of the moments of the order parameter 
do not show the characteristic signatures that should be present if of a low-temperature 
reentrant transition occurred. Our results suggest that renormalized spinwave theory is 
applicable as T — > for sufficiently weak disorder. This is inconsistent with theories that 
predict the reentrance of the paramagnetic phase at low temperature JTJ but agrees with 
more recent theories [IU|,n . 
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FIGURES 

FIG. 1. Schematic phase diagram of the XY model with random phases in the T — a plane 
according to the theory of Rubinsteinei al. T + and T_ are, respectively, the values of the KT 
transition temperature and of the reentrance temperature for the case of an infinite vortex core 
energy. 

FIG. 2. The Binder function as a function of temperature for different system sizes and several 
values of the strength of the disorder in the weak disorder regime, (a) a = 0.09, (b) a = 0.16, (c) 
a = 0.25, and (d) a = 0.36. 

FIG. 3. A scaling plot of the data presented in the previous figure. For each value of a the 
data for all system-sizes and temperatures collapse into a single function of L/£(T,a). 

FIG. 4. The Kosterlitz-Thouless transition temperature as a function of disorder as determined 
from the scaling plots. The dotted line is the prediction of the modified RG equations. The 
temperatures at which reentrance was expected to occur are indicated by the open circles. 

FIG. 5. The fourth moment of the order-parameter as a function of temperature for a = 0.25 
and several system sizes. 

FIG. 6. Logarithmic plot of the fourth moment of the order-parameter as a function of sys- 
tem-size for o = 0.25 and several temperatures. The latter go between 0.1J and 0.7J in steps of 
0.05J. The top curve corresponds to the lowest temperature. 

FIG. 7. The correlation function exponent as a function of temperature for several values of a. 

FIG. 8. The Binder function as a function of temperature for different system sizes and several 
values of the strength of the disorder in the strong disorder regime, (a) a = 0.49, (b) a = 0.64, (c) 
a = 0.81, and (d) a = 1. The solid lines are fits to the functional form described in the text. 



17 



FIG. 9. The Binder function as a function of temperature and disorder for several values of the 
strength of the disorder and L = 16. 



FIG. 10. The universal scaling function in the disordered phase as a function of x = L/£. 

FIG. 11. Temperature dependence of the correlation length (in arbitrary units) in the strong 
disorder regime for the same values of a as in Fig. 9. The highest curve corresponds to the lowest 
value of the disorder. The full lines are derived from the fits in Fig. 8. 

FIG. 12. The zero-temperature correlation length as a function of disorder. The unit of length 
is arbitrary. The solid line is a fit to the expression £ ~ exp[6/(cr c — a)] with b = 0.69. The 
divergence at a c signals the disorder-driven unbinding of vortices. 
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